/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "emptyPolyPatch.H"
#include "simpleControl.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"

    Info << "Reading field p" << nl << endl;
    volScalarField p
    (
        IOobject
        (
            "p",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );
    
    forAll(mesh.boundaryMesh(), patchI)
    {
        const polyPatch& meshPatch = mesh.boundaryMesh()[patchI];
        word name = mesh.boundaryMesh()[patchI].name();
        if(!isA<emptyPolyPatch>(meshPatch) && !mesh.magSf().boundaryField()[patchI].coupled()) 
        {
            forAll(meshPatch, pfaceI)
            {
                doubleScalar xx = M_PI * mesh.boundary()[patchI].Cf()[pfaceI].x();
                doubleScalar yy = M_PI * mesh.boundary()[patchI].Cf()[pfaceI].y();
                p.boundaryFieldRef()[patchI][pfaceI] = Foam::cos(xx) * Foam::cos(yy);
            }
        }
    }

    volScalarField p_e
    (
        IOobject
        (
            "p_e",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        p
    );
    
    Info << "Reading transportProperties\n" << endl;
    IOdictionary transportProperties
    (
        IOobject
        (
            "transportProperties",
            runTime.constant(),
            mesh,
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    );

    Info << "Reading the diffusivity\n" << endl;
    dimensionedScalar gamma (transportProperties.lookup("gamma"));

    volScalarField omega(p);
    forAll(omega,cellI)
    {
        doubleScalar xx = M_PI * mesh.C()[cellI].x();
        doubleScalar yy = M_PI * mesh.C()[cellI].y();
        omega[cellI] = -2 * M_PI * M_PI * Foam::cos(xx) * Foam::cos(yy);
    }
    simpleControl simple(mesh);
    // while (simple.correctNonOrthogonal())
    // {
    //     fvScalarMatrix pEqn
    //     (
    //         fvm::laplacian(gamma, p) == omega
    //     );

    //     label pRefCell = 0;    
    //     doubleScalar pRefValue = 0.0;
    //     pRefCell = mesh.findCell(vector(0.5,0.5,0.05));
    //     pEqn.setReference(pRefCell,pRefValue);

    //     pEqn.solve();
    // }
    int i = 0;
    while (i<50)
    {
        fvScalarMatrix pEqn
        (
            fvm::laplacian(gamma, p) == omega
        );

        label pRefCell = 0;    
        doubleScalar pRefValue = 0.0;
        pRefCell = mesh.findCell(vector(0.5,0.5,0.05));
        pEqn.setReference(pRefCell,pRefValue);

        pEqn.solve();
        forAll(p,cellI)
        {
            p[cellI] = 0.0;
        }
        i++;
    }

    Foam::Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    
    doubleScalar L1 = 0.0, Lu1 = 0.0, Ll1 = 0.0;
    doubleScalar L2 = 0.0, Lu2 = 0.0, Ll2 = 0.0;
    forAll(p_e,cellI)
    {
        doubleScalar xx = M_PI * mesh.C()[cellI].x();
        doubleScalar yy = M_PI * mesh.C()[cellI].y();
        p_e[cellI] = Foam::cos(xx) * Foam::cos(yy);
    }
    forAll(p,cellI)
    {
        Lu1 += mag(p[cellI] - p_e[cellI]);
        Ll1 += mag(p_e[cellI]);
        Lu2 += sqr(p[cellI] - p_e[cellI]);
        Ll2 += sqr(p_e[cellI]);
    }

    reduce(Lu1,sumOp<doubleScalar>());
    reduce(Ll1,sumOp<doubleScalar>());
    reduce(Lu2,sumOp<doubleScalar>());
    reduce(Ll2,sumOp<doubleScalar>());

    L1 = Lu1 / Ll1;
    L2 = Foam::sqrt(Lu2 / Ll2);
    Info << "L1: " << L1 << nl;
    Info << "L2: " << L2 << nl;

    std::ofstream outfile;
    outfile.open("error.txt", std::ios::out | std::ios::app);
    outfile << mesh.C().size() << "\t" << L1 << "\t" << L2 << std::endl;
    outfile.close();

    p.write();
    p_e.write();
    Info << nl << "End" << nl << endl;
    return 0;
}

